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Abstract. We present an event-by-event study of cosmic ray (CR) composition with the 
reflected Cherenkov light method. The fraction of CR light component above 5 PeV was 
reconstructed using the 2013 run data of the SPHERE experiment which observed optical 
Vavilov-Cherenkov radiation of extensive air showers, reflected from snow surface of Lake Baikal. 
Additionally, we discuss a possibility to improve the elemental groups separability by means of 
multidimensional criteria. 


1. Introduction 

The study of the superhigh energy {E > lO^^eV = 1 PeV) cosmic ray (CR) composition is the 
most difficult experimental problem of CR physics. Indeed, even for the most basic primary 
nuclei mass-dependent quantity, the mean logarithmic mass number < In A >, the spread of 
results is very high (almost from proton to Iron at some energies) [1]. A somewhat newer 
review, at the cost of ignoring a part of older works, shows lesser, but still considerable, scatter 
of < In A > values [ 2 ]. Available information on the elemental groups spectra [ 3 ]-[ 7 ] at £’ > 3 
PeV is almost exclusively due to particle detectors relying on the electron and muon numbers 
measurement [ 2 ]; this method is highly model dependent [H]. Moreover, results obtained with 
event-by-event techniques are very scarce [2], and most of the above-mentioned works 0 - 11 ! 
relied on the deconvolution method, that may suffer from the ill-posedness of the inverse problem 

m- 

In the present work we describe an event-by-event study of CR composition using the 2013 
run data of the SPHERE experiment. The results of this Cherenkov experiment are much 
less dependent on high-energy hadronic model than for the case of particle detectors and, 
additionally, are expected to be more stable than ones obtained with the deconvolution method. 
The datasample used in our analysis is considered in section 2. The method employed for the 


composition study is described in section 3. In section 4 the result of the analysis — the fraction 
of CR light component above 5 PeV — is presented. Finally, in section 5 we show that it 
is possible to enhance the sensitivity of the method to the primary nuclei mass number with 
multidimensional criteria. 

2. The datasample and low-level data analysis 

The SPHERE experiment, stage 2 (2008-2013), employed the SPHERE-2 balloon-borne detector 
HU to observe optical Vavilov-Cherenkov radiation (’’Cherenkov light”) of extensive air showers 
(EAS), reflected from snow surface of Lake Baikal. The foundations of the method were 
discussed in m- Many details on the detector hardware, observation conditions, simulations, 
and experimental data analysis could be found in [TT]; for brief review see [13]. During the 2013 
run, the typical observation altitude (measured by GPS) was in range H— 400-700 m above the 
snow surface. In total, 3813 events were recorded; 459 of them were recognized as EAS events 
(for the low-level experimental data analysis methods, see [H], section 5.4). For the next step 
of the analysis, the lateral distribution function (LDF) reconstruction, we employed a new code, 
that is generally similar to the code used in our previous works HU, HU, but was written in 
a completely independent manner. This circumstance allows us to test the robustness of our 
results against those obtained with the ’’old” LDF reconstruction method HU- 

In total, 421 LDF were obtained with the new method out of the 459 events. The zenith and 
azimuthal angles of these showers {O^ip) were measured under an assumption that the shower’s 
front is a plane. Out of the 421 LDF 354 events have estimated zenith angles 9 < 40° and 
reconstructed axis distance to the center of the detector’s fleld-of-view (FOV) projection to 
the observation surface R < RMax] RMax = Rfov+30 m, where Rfov = 0.4903 • H is the 
FOV projection radius. The present analysis relies on a subsample of the late 354 LDF, 328 
events. The simulation, however, was perfomed for the full sample of showers with 9 < 40° 
and R < Rjviax^ that corresponds to 354 events. Given that the experimental sample is nearly 
complete, in what follows we neglect all possible associated biases. 

3. The method to separate elemental groups 

3.1. The LDF steepness parameter 

For composition study we follow the general approach developed in HU-HU, HU- Using a 
sample of model LDF, we deflne the criteria of the elemental groups separation for different 
primary energies, zenith angles and observation altitudes. As in HU-HU, HU, the parameter 
sensitive to the primary nuclei mass is the LDF ’’steepness”, i.e. the ratio of the number of 
Cherenkov photons in the circle with the radius of 70 m with the center in the LDF’s axis to 
the same number in the concentric ring with the radii of 70 m and 140 m. 

3.2. Simulations 

Simulations of EAS development were carried out for primary energies Emc= 10, 30, 100 PeV 
by means of full direct Monte Carlo (MG) method using the CORSIKA code [T7] with the 
QGSJET-I high energy hadronic model m and the GHEISHA low energy hadronic model HU- 
The detector response database, that consists of a large number ('^ 3-10^) of model showers, was 
calculated using the Geant4 code m- Each ’’CORSIKA shower” was used multiple times with 
different axis coordinates. The axis coordinates were uniformly (randomly) distributed over big 
square with dimensions 1.5 • H x 1.5 • H. 

As well, the instrumental trigger response was simulated for a range of energies Ethq =5-200 
PeU, and for a wide range of other conditions. For Ethq =5-17.3 PeU, the 10 PeU CORSIKA 
showers were used, and for Ethq >17.3 PeU — the 30 PeV showers. The model showers with 
Errig/E mc = K were obtained by multipying the corresponding responses to the factor K. We 



Figure 2. Reconstructed energy for model 
showers with various primary nuclei, primary 
energies and observation altitudes. 


Figure 1. An example of the individual 
model LDF (red circles) together with a 
’’composite model LDF” (green circles) that 
fits the individual LDF. The energy of 
primary nucleus was set to 100 PeV. The 
legend shows the values of some ’’true MC” 
parameters; (xo,2/o) ['^] is the axis position 
coordinates with respect to the center of the 
detector’s FOV projection to the observation 
surface. The zenith angle value 9 is measured 
in °. As well, the reconstructed energy value 
(96.18 PeV) is shown. 



Figure 3. The same, as in figure 2, but only 
for showers that were registered by the model 
of the instrumental trigger. 


have verified that the Emc values discreteness has a negligible impact on results of the trigger 
response simulation. 

3.3. Estimation of shower parameters 

As a first step of the elemental groups separation procedure, we have estimated the LDF 
steepness parameter r] for all model showers, integrated on time, by fitting these ’’discrete LDF” 
by a set of specially calculated smooth continious curves. Each curve of the set (we call it a 
’’composite model LDF”, CLDF) was composed from many model ’’discrete LDF” corresponding 
to the same CORSIKA shower. An example of a CLDF is shown in figure 1 by green circles, and 
the corresponding individual model LDF is plotted by red circles. The normalisation factor of a 
CLDF, that is nearly proportional to the primary energy EH, was utilized as the primary energy 
estimator, and the CLDF steepness parameter — as an estimate of the r] value. The regions (0, 
70) m and (70, 140) m used for the r] parameter estimation are constrained by vertical black 
lines in figure 1. 

The energy reconstruction performances of the described procedure are presented in figures 
2-3. Figure 2 deals with all showers in the model sample, irrespectively of the axis position. The 
mean reconstructed energy values are shown by symbols for proton and Iron primaries for the 
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Figure 4. The result of the model showers 
classification for Ethq— 11*2 PeV. 
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Figure 5. The same, as in figure 4, but for 
Exrig^ 15.9 PeV. 
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Figure 6. The same, as in figures 4-5, but 
for ETrig= 20.9 PeV. 


Figure 7. The same, as in figures 4-6, but 
for ETrig= 29.6 PeV. 
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Figure 8. The same, as in figures 4-7, but 
for Errig^ 41.8 PeV. 
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Figure 9. The same, as in figures 4-8, but 
for ETrig— 102.6 PeV. 






































primary (’’true MC”) energy values Emc— (10, 30, 100) PeV^ and several altitudes (400, 500, 
580, 700, 900) m. The corresponding standard deviations ge {cte is a measure of the primary 
energy statistical reconstruction uncertainty) are shown by bars. A small artificial shift between 
the results for different conditions is introduced on the H axis to make them visible. Similar 
results for Emc= 10 PeV and 30 PeV are shown in figure 3 for showers that satisfy the trigger 
condition. We do not show results for 900 m, because the maximal altitude for the 2013 run 
was ?=^700 m. It is evident that the statistical uncertainty of the primary energy value is qiute 
low (from ^10 % to ^20 % depending from the primary energy and the altitude values), and 
the systematic error is moderate (<5 % for most of the cases). 

3.4‘ Separation of the elemental groups 

To study the CR composition, we employed a set of Bayesian classifiers [22] using the above- 
described model LDF sample as a training set. In the present work, for simplicity, we constrain 
ourselves to the case of the linear border between the classes of primary nuclei, and two event 
classes — proton and Iron primary nuclei. Following the famous ’’Bayes rule”, we have set the 
prior probabilities 0.5 for both classes. The results of the classifier training are shown for the 
altitude value H= 400 m and different energies Ethq in figures 4-9. The values of the LDF 
steepness parameter for proton (red crosses) and Iron (blue circles) vs. the primary zenith angle 
value 9 are plotted; the boundary between the classes is shown by green straight line. The Emc 
parameter value is shown (see subsection 3.1 for explanation), as well as the parameters of the 
border line equation. The numbers of proton (Np^up and Np) and Iron {Nf_up and Nf) showers 
above the border, and total numbers of these model showers, respectively, are presented. To 
suppress statistical fluctuations arising from a small number of CORSIKA showers, we have 
introduced a small Gaussian smearing on both parameters (^, 77 ). 

The fraction of model proton showers classified as light nuclei is Ep = Np^up/Ep— (0.635, 
0.519, 0.374, 0.359, 0.335, 0.350) for the Erng values ( 11 . 2 , 15.9, 20.9, 29.6, 41.8, 102.6) PeV, 
respectively, and the corresponding fraction of Iron contamination Ef — Nf_y^p/Nf— (0.050, 
0.025, 0.007, 0.008, 0.007, 0.009). It is evident from figures 4-5 that Iron nuclei are effectively 
suppressed by the instrumental trigger effects; at sufficiently low energy the SPHERE-2 detector 
could observe practically only light nuclei. Overall, the Iron contamination appears to be quite 
low. Results such as shown in figures 4-9 were obtained for 5 altitudes (400, 500, 580, 700, 900) 
m and 54 logarithmically spaced energy bins from 5 to 200 PeV. 

All experimental showers from the sample defined in section 2 have the parameters {9,E,r]) 
estimated, and the altitude value H measured. These events were classified as proton-like or 
Iron-like according to their position above or below the border. The total number of the light 
CR component showers in 6 energy bins in the energy range 5-200 PeV was reconstructed as 
Efi-iight = Efi_p_iikelEp, where Ni_p_iike is the number of experimental showers above the border 
in the ith bin. The light and heavy component intensities were corrected for the instrumental 
acceptance effects to account for the difference between the detection efficiencies of light and 
heavy nuclei. 

In comparison to |14] , in the present work we were able to substantially lower the threshold of 
the method using calculations of the instrumental acceptance mm- The first-order correction 
in the present work uses the 2013-integrated acceptance curves for the case of i? < Reov + 100 
m, instead of i? < Reov + 30 m. The number of available model showers at E— 100 PeV is 
currently low, therefore the analysis was performed with the 10 PeV and 30 PeV model showers. 
Corrections for these inaccuracies are made in section 4. 

4. Results 

The result of our analysis — the estimated fraction of light nuclei in the energy range 5-200 
PeV — is shown in figure 10 by red stars with statistical uncertainties. For the case of the 




Figure 10. The fraction of light (’’proton¬ 
like”) nuclei reconstructed using the 2013 run 
data of the SPHERE experiment. 


Figure 11. The fraction of light nuclei before 
(red stars) and after corrections (black and 
blue circles; see text for more details). 


120-240 PeV energy interval the statistical uncertainty is quite high and our results are not 
constraining. For the case of 3.5-7 PeV energy interval we didn’t find any indication of the Iron 
nuclei presence. To estimate the statistical uncertainty for this bin, we have added one ’’bonus” 
Iron shower; the corresponding result is shown in figure 10 by the blue star. The light nuclei 
fraction behaviour presented in figure 10 is qualitatively similar to the KASCADE-Grande result 
[7], and doesn’t contradict our earlier estimate [H] obtained with the 2012 run data. We have 
already applied the correction for model energy discreteness here (see section 3). 

Several systematic effects do exist that could change our results. A shift of the light nuclei 
fraction could occur as a result of the energy estimation systematic uncertainty (see figure 3). 
For some energy bins, this effect tends to enhance the reconstructed light nuclei intensity up to 
10-15 % with respect to its true value. Therefore, the corrected fraction of light nuclei in these 
bins would be lower than plotted in figure 10 (for instance, 25-27 % instead of 30 %). We have 
estimated and applied corrections for this effect, as well as the acceptance correction discussed 
at the end of section 3. The new result after these corrections is shown in figure 11 by black 
circles. As well, we have estimated the LDF reconstruction procedure systematic uncertainty, 
and the corresponding result including this last correction, as well as all previous ones, is shown 
in figure 11 by blue circles. Dedicated work to study these and other systematic effects is in 
progress, and will be reported elsewhere. 

5. Multidimensional method for superhigh energy CR composition study 

In sections 3-4 were have employed the set of one-dimensional criteria, relying on the r] = 
r]{A^ H) parameter. Here we show that it is possible to improve the nuclei classes 

separability by using multidimensional pattern recognition methods. We demonstrate the 
case on our model sample (without imposing the trigger effects) by studiyng classification 
performances for 3 nuclei groups separation — proton. Nitrogen, and Iron. The feature vector 
was chosen to be 4-dimensional, and was composed of Cherenkov light intensities in concentric 
rings with the center in the shower’s axis and with radii (0, 40) m, (40, 80) m, (80, 120) m, 
and (120, 160) m. The zenith angle value 9 was added as the fifth classification parameter. 
The multivariate Bayesian pattern recognition technique under assumption of multidimensional 
Gaussian distribution of features m was used for this study. This technique was recently shown 
to be very effective for the case of CR composition study with lateral-angular distribution of 









Cherenkov light m-m- 


5.1. Proton-Nitrogen separation 

All graphs in this subsection (figures 12-15) show contamination of Nitrogen nuclei (i.e. the 
fraction of Nitrogen nuclei classified as protons) vs. proton selection efficiency. Black curves are 
drawn for H— 400 m, red curves — for H— 580 m, and blue curves — for H— 900 m. The 
primary energy values and zenith angle range are shown in the captions to the graphs. 




Figure 12. Proton-Nitrogen separation for Figure 13. The same, as in figure 12, but for 
10 PeV and (9= 0-20°. 6>= 20-40°. 




Figure 14. The same, as in figures 12-13, 
but for 30 PeV and 9= 0-20°. 


Figure 15. The same, as in figure 14, but for 
(9= 20-40°. 


5.2. Proton-Iron, proton-Helium and Nitrogen-Iron separation 

The same results as in figures 12-15, but for separation of other elemental groups, were also 
obtained. For the case of proton-iron separation (except the 10 PeV and 0= 20-40° case), 
the typical proton selection efficiency corresponding to ^1 % Iron contamination is ?^60-90 %, 
depending on conditions. For the 10 PeV and 9= 20-40° case, the proton selection efficiency 
is lower, but in this case the Iron showers would be suppressed by the trigger effects, as was 
explained in section 3. The same applies to other cases where proton selection efficiency is <75 
%; therefore, the proton selection efficiency is, in fact, >75 % for ^1 % Iron contamination. 

As well, we have studied the Nitrogen-Iron classification performance, and obtained the 
typical Nitrogen selection efficiency value ~40 % for ^1 % Iron contamination. To conclude. 






let us note that proton-Helium separation is, in principle, also possible with the presented 
multidimensional technique. However, with the criteria considered in the present work, the 
proton selection efficiency was observed to be always <25 % for ^1 % Helium contamination. 

6. Conclusions 

We have presented and discussed the method for event-by-event study of superhigh CR 
composition applicable to the SPHERE experiment working conditions. It was found that 
a part of the light component showers can be readily selected against the heavy component 
background. The primary energy can be reconstructed with sufficiently low uncertainty, both 
statistical and systematic. The result of the analysis — the light nuclei fraction vs. energy — 
is qualitatively similar to the KASCADE-Grande result, and doesn’t contradict to our earlier 
estimates. As well, we have shown that multidimensional criteria, in principle, allow to enhance 
the separability of the elemental group classes. 
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